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We  consider  Monte  Carlo  simulation  of  transient 
behavior 'of  continuous^ time  Markov  processes  with  large 
finite  state  spaces.  The  usual  event-scheduling  ap¬ 
proach  is  compared  to  four  alternative  approaches 
which  do  not  use  event  lists.  Variance  reduction, 
programming  effort,  and  storage  requirements  are  con¬ 
sidered.  All  Monte  Carlo  simulation  programs  are 
written  in  SIMSCRIPT  II. 5.  These  approaches  are  then 
compared  to  a  deterministic  method  called  Randomiza- 


The  four  alternative  Monte  Carlo  simulation  meth¬ 
ods  are  based  on  two  dichotomies.  The  first  concerns 
continuous  versus  discrete  time:  the  process  can  be 
simulated  directly  in  continuous  time  or  the  uniform- 
ized  embedded  discrete^t ime  Markov  chain  can  be  simu¬ 
lated  and  then  analytically  randomized  to  recover  the 
original  continuousf time  process.  The  second  dichoto¬ 
my  concerns  a  table-flook?up  versus  an  algorithmic  ap¬ 
proach:  the  mean  holding  time  and  the  transition  dis¬ 
tribution  for  the  current  state  of  the  simulated  pro¬ 
cess  can  be  read  from  memory  (requiring  large  memory 
storage)  or  computed  (requiring  more  computation  time). 

The  above  methods  are  applied  to  the  simulation  \ 
of  a  multi-echelon  repairable  item  provisionong  system  1 
with  a  finite  population  of  repairable  items  and  lim¬ 
ited  repair  capacity.  Various  systems  are  used.  Prob¬ 
abilities  of  system  availability  are  estimated.  The 
variances  of  these  estimates  are  multiplied  by  CPU 
times  to  get  measures  of  computational  efficiency  of 
the  different  simulation  methods. 
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The  q..'s  are  the  transition  rates.  The  infinitesimal 
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generator  Q  seems  to  be  the  most  natural  way  to  de¬ 
scribe  the  stochastic  nature  of  continuous  time  Markov 
models  with  denumerable  state  spaces.  The  state  prob¬ 
ability  vector  at  time  t  is  denoted  jr ( t )  =  £tt  (t), 

7T2(t) ,  ...»  TT^(t)),  where  tt  (t)  =  P^X(t)  =  s},  s  E  S . 

These  transient  probabilities  satisfy  the  Kolmogorov 
forward  equation 

JT’(t)  =  TT(t)Q,  t  >  0.  (1) 

This  is  an  initial  value  system  with  £(0)  given.  (See 
[1,2,3]  for  more  background  on  Markov  processes.) 

The  purpose  of  this  paper  is  to  consider  differ¬ 
ent  ways  of  computing  or  estimating  the  above  tran¬ 
sient  probabilities,  £(t),  t  >  0,  i.e.,  solution  of 
the  Kolmogorov  equation  (1).  There  are  two  general 
approaches:  deterministic  algorithms  and  Monte  Carlo 
algorithms.  We  shall  consider  one  deterministic  algo¬ 
rithm:  the  Randomization  Technique.  We  shall  consid¬ 

er  five  Monte  Carlo  approaches  for  simulating  the 
transient  behavior  of  continuous  time  Markov  processes. 

THE  RANDOMIZATION  ALGORITHM 


INTRODUCTION 

Let  { X C t ) ,  t  ^  0}  be  a  continuous-time  time- 
homogeneous  Markov  process  on  a  finite  state  space  S  = 
{l,2,...,m}.  All  such  Markov  processes  can  be  charac¬ 
terized  by  an  initial  distribution  J.(0)  and  an  infini¬ 
tesimal  generator 
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Any  Markov  process  X  on  a  finite  state  space  can 
be  represented  as  a  discrete  time  Markov  chain  (the 
uniformized  embedded  chain)  "randomized"  by  a  Poisson 
process.  Define 

P  =  Q/A  +  I,  where  A  =  max  q.  (2) 

icS  1 

and  1  is  the  identity  matrix;  P  is  a  stochastic  matrix. 

Let  {Y  ,  n  =  0,1,2,...}  be  a  Markov  chain  on  S  with 
n 

transition  matrix  P  and  initial  distribution  ^(0) . 

Let  { N ( t ) ,  t  ^  0l  be  a  Poisson  process  with  rate 
which  is  independent  of  { ,  n  -  0,1,2,... Then 

^ Y n ( t ) •  t  ?  0  is  a  Markov  process  with  generator  Q 

and  initial  distribution  _^ ( 0 )  and  hence  is  probabilis¬ 
tically  identical  to  ;X(t),  t  >  0‘.  This  construction 
makes  it  possible  to  compute  transient  probabilities 
of  a  Markov  process  with  generator  Q  from  transient 
probabilities  of  a  Markov  chain  Y  with  transition  ma¬ 
trix  P  and  a  Poisson  process  N  with  rate  The  tran¬ 

sient  prohabi  1  i  t  ies  ,>t  V  are  denoted  ;(n)  = 

1 

!  ,  ( n ) ,  ....  ( n  0  .  whe re  .*  (  n  >  =  P  ( Y  -  s  >  ,  s  •  5 .  The 
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or  equivalently, 
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See  Cross  and  Miller  [4J  for  additional  discussion  and 
details.  (Equation  (3)  can  also  be  found  in  Qinlar 
[L,  p.  259].) 


The  infinite  series  in  Equation  (3)  must  be  trun¬ 
cated  for  computational  purposes.  Let 


T(e,t)  =  min 


k: 
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where  c  equals  an  acceptable  error  (specified  by  the 
user).  The  computational  version  of  Equation  (3)  is 


i  (t) 


T(e,t) 

I  £(n) 

n=0 


-At 


(At)11 


(5) 


Truncation  of  the  infinite  series  involves  a  probabil¬ 
ity  loss  of  at  most  £;  thus  all  probabilities  (of 
states  or  subsets  of  states)  will  have  an  error  be¬ 
tween  -e.  and  0.  Note  that  the  randomization  formula 
(5)  reduces  the  calculation  of  transient  probabilities 
of  a  Markov  process  to  those  of  a  Markov  chain  and  un¬ 
derlying  Poisson  process,  both  of  which  are  more  amen¬ 
able  to  exact  numerical  evaluation. 

The  ti>’s  are  computed  recursively  using  the  rela¬ 
tion  from  standard  Markov  chain  theory: 

<fc(0)  =  *.(0);  £(n+l)  =  £(n)P,  n  >  0.  (6) 

(Note  that  Equation  (6)  involves  only  nonnegative  num¬ 
bers,  a  fact  that  contributes  to  numerical  stability 
of  the  algorithm.)  The  matrix  P  is  usually  sparse  and 
thus  the  above  matrix  multiplication  should  be  per¬ 
formed  by  an  appropriate  algorithm.  Such  a  multipli¬ 
cation  algorithm  is  described  by  Cross  and  Miller  [41. 
The  number  of  operations  in  this  algorithm  is  propor¬ 
tional  to  the  sum  of  the  number  of  states  and  the  num¬ 
ber  of  transitions. 


In  short,  the  standard  randomization  computation¬ 
al  algorithm  computes  A  and  P  from  the  generator  Q  us¬ 
ing  (2).  It  computes  the  truncation  point  T(i:,t)  from 
(4),  then  the  ^>(n)'s  using  (6)  recursively,  accumulat¬ 
ing  in  Equation  (5)  to  give  jT‘(t). 


Cross  and  Miller  [5]  have  computed  transient  prob¬ 
abilities  for  Markov  process  models  of  multi-echelon 
inventory  systems  with  20,000  states  and  200,000  tran¬ 
sitions  using  the  randomization  algorithm.  Melamed 
and  Yadin  (61  have  applied  the  method  to  queuing  net¬ 
works  with  a  large  number  of  states.  Miller  (7)  has 
adapted  the  randomization  algorithm  to  efficiently 
handle  certain  stiff  systems  which  arise  in  the  reli¬ 
ability  analysts  of  fault-tolerant  systems. 


There  are  other  deterministic  approaches  to  the 
solution  of  the  Kolmogorov  equation  (1)  which  we  will 
not  consider.  Two  general  approaches  are:  (I)  numer- 
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leal  integration  techniques  such  as  Runge-Kutta, 
predictor-corrector,  etc.;  and  (ii)  exponentiation 
Qt 

[ tt ( t )  =  jr(0)e  ]  by  computing  the  spectrum,  computing 
the  Taylor  series,  or  other  means.  The  randomization 
technique  has  a  distinct  advantage  over  these  ap¬ 
proaches  in  that  a  bound  on  the  global  error  can  be 
set  by  the  user,  and  it  is  achieved.  (The  only  other 
source  of  error  is  the  influence  of  rounding  and  trun¬ 
cations  by  the  machine  performing  the  calculations;  by 
noting  that  the  randomization  algorithm  mainly  in¬ 
volves  multiplication  and  addition  of  positive  num¬ 
bers,  Grassnann  [8]  has  bounded  this  error.)  Further¬ 
more,  Grassmann  [9J  has  shown  empirically  that  random¬ 
ization  is  more  efficient  for  some  queuing  systems. 

We  use  the  randomization  algorithm  as  a  determin¬ 
istic  method  for  solving  the  Kolmogorov  equation  (1). 
We  also  use  the  concept  of  randomizing  a  discrete-time 
Markov  chain  (as  described  above)  as  part  of  a  Monte 
Carlo  simulation  method. 


MONTE  CARLO  SIMULATION  METHODS 

We  use  five  Monte  Carlo  simulation  methods.  The 
first  uses  the  classical  world  view  of  "future  event 
scheduling."  The  others  exploit  the  Markovian  nature 
of  the  system  and  use  no  future  event  list. 

Event-scheduling  Method  (ES) 

The  classical  event-scheduling  method  is  well 
known.  It  is  used  by  languages  such  as  SIMSCKIPT.  An 
event  is  defined  as  a  change  of  state  of  the  system. 

The  occurrence  of  one  event  often  triggers  the  occur¬ 
rence  of  another  event  in  the  future;  thus,  when  the 
first  event  occurs  the  simulation  algorithm  schedules 
the  future  occurrence  of  the  other  event  by  creating 
an  event  notice  and  filing  it  in  a  "future  event  list." 
Advancing  simulated  time  requires  searching  the  event 
list  for  the  next  scheduled  event  and  noting  its  time 
of  occurrence.  The  event-scheduling  approach  is  very 
general;  it  is  applicable  to  virtually  all  discrete 
event  systems. 

Simulating  Markov  Processes 

If  the  system  is  Markovian,  there  may  be  advan¬ 
tages  in  exploiting  that  property.  An  event-scheduling 
algorithm  can  be  used,  but  it  may  consume  considerable 
execution  time  for  filing  and  retrieving  notices  from 
the  future  event  list  if  the  simulated  system  is  even 
slightly  complex.  It  may  be  more  economical  to  simu¬ 
late  a  Markovian  system  directly  as  a  Markov  process. 
There  are  various  wavs  to  do  this,  for  example,  see 
Hordijk,  Iglehart,  and  Schassberger  [10]  or  Fishman 
nil. 


Consider  a  continuous  time  Markov  process 
{X(t),  t  >  0}  with  initial  probability  vector  _(0)  and 
infinitesimal  generator  Q.  Let  equal  the  time  of 

the  n-th  transition  of  X,  n  =  1,2,...,  I’  =  0.  The 

process  [7  -  X(T  ),  n  =  0,1,2,...-  is  an  embedded 

n  n 

discrete-tine  Markov  chain  (jump  chain).  The  transi¬ 
tion  matrix  of  7.  is 
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where  =  q ^ j /q ,  i  j*  j  ,  The  continuous- time  Ilarkov 

process  X  is  really  a  semi-Markov  process  (see  [1,2, 

3))  with  transition  matrix  R  and  exponential  holding 

times:  P(T  -  T  >  t  I  Z  =  i)  =  exp(-q.t);  i  a  L, 
n+1  n  '  n  r  T  ’ 

n  =  0,1,2,....  This  structure  will  be  used  in 
two  Monte  Carlo  simulation  algorithms. 

Table-look-up  Continuous- tine  Method  (TC) .  In 
this  Monte  Carlo  algorithm  we  use  the  structure  of  em¬ 
bedded  (jump)  chain  and  exponential  holding  times. 

The  transition  matrix  R  of  the  jump  chain  and  the  vec¬ 
tor  of  rates  (q j »^2* *  *  * »^n)  are  stored.  To  simulate 

one  sample  path  {x(t),  0  <  t  <  si  of  X  we  proceed  as 
fo 1 lows : 

(i)  let  t  =  0; 

(i£)  compute  x(t),  a  random  variate  from  the  prob¬ 
ability  vector  jl(0) ,  using  the  inverse  cdf  meth¬ 
od  ,  see  (til; 

(iii)  generate  j,  a  random  variate  whose  probability 
vector  is  the  x(t)-th  row  of  R,  using  the  in¬ 
verse  cdf  method; 

(iv)  generate  an  exponential  (rate  q  /t\)  random  var¬ 
iate  and  add  it  to  t; 

(v)  let  x(t)  ■  j; 

(vi)  repeat  steps  (iii)  through  (v)  until  t  >  s. 

This  TC  algorithm  has  the  advantage  of  being  quite 
simple  and  fast.  It  has  the  disadvantage  of  requiring 
storage  of  all  nonzero  elements  of  R  and  the  vector  of 
rates,  =  (q ? , . . . ,qn) •  Also,  complex  systems  may 

have  generators  which  require  sophisticated  programs 
for  their  computation. 

Algor i thnic  Cont inuous-t ime  Me thod  (AC) .  This 
algorithm  is  the  sane  as  TC  except  instead  of  storing 
the  rate  vector  £  and  all  the  nonzero  transition  prob¬ 
abilities  in  R,  an  algorithm  computes  the  values  as 
they  are  needed.  In  this  approach  the  state  x  of  the 
system  (which  itself  may  be  a  vector)  Is  more  descrip¬ 
tive  than  with  TC,  where  states  are  simply  numbered; 
the  state  descriptions  x  lor  AC  are  the  same  as  or  sim¬ 
ilar  to  KS.  Suppose  there  are  at  most  k  events  that 


KV ( ,  KV 


To  simulate  a  sample 


pa tli  using  AC  we  proceed  as  follows: 

(i)  let  i  =  0; 

(ii)  generate  an  initial  state  and  let  x(t)  3  x^; 

(iii)  compute  rales  ot  occurrence  q^x(t),ij  of  event 
KVj,  i  =  I.  2,...,  k,  given  the  current  state 

x  ( t  )  ,  let 


j(x(t  ) ,  ij,  i  i  xtt  ) ; 


generate  ti>e  next  event  to  i<>  >  ur  using  a  Monte 
Carl«»  onpiil  at  i  on  on  the  state  des,  ript*>r  •:  ( t  »  , 
then  *  . » 1  I  the  appropriate  event  routine  (simi¬ 
lar  to  IS)  to  <  ampul'-  the  ■  ll.lll’r  'I  sf  i  (  ,  -  dr*»- 
*  f  1  J  •  t  •  ■  I  l  • '  i'et  a  UeW  >t  it  e  •.  ; 

.Tiicf  i f  < •  n  f-.t u-neiti  i  i.'  (t  it  i  u  )  t  i a,/. " . 

•  •  (  t  i 

at  i  it «  nu!  »d.i  it  t .  •  t  ; 


I  1  I  ill' t  i  I  t 


At  f  i  a  s  t  I  i  >  r,:.  11  f  i'<  '  t  «  ;i.  *  J  j  v. 

t  i f  l  a i  it  i -  1 1  I s  t  r  i n ■ .  i  i  j  • .  *  .  \  ; •  i  .i 

ha-,  tin  i<!\  iM  i-t  t  a-.i,  i  :■  •  -e-.  -  f  .  ■  . j  >  1 1  r . 
-apahilit.  t  hand;  lu  mil",  it*  .1,1.  .p 
imp  l*-i  1 1  m  - 1 1  i  :  .  .i  vm  •  i  i  - 
I'hvs  ).  I  ;  ■  1  'pi  r  !  l  I  .  !  t  \|  .  !  I  '  : . 


Table-look-up  Discrete-time  Method _(TD) .  This 
method  is  the  sane  as  TC  except  we  simulate  the  utii- 

fornized  embedded  discrete- time  Markov  chain  {Y  , 

n 

n  =  0,1,...)  which  has  initial  distribution  j(0)  and 
transition  matrix  P  =  Q/A  +  f;  see  Equation  (2).  We 
then  convert  this  to  the  continuous-time  process 
(X(t),  t  ^  0}  by  randomizing,  i.e.,  using  Equation  (5). 
In  particular,  we  simulate  the  sample  path  {y  ,  n  *  0, 
l,2,...,n}  as  follows:  1 

(i)  let  n  *  0; 

(ii)  compute  yn,  a  random  variate  from  the  probabil¬ 
ity  vector  jr(0)  ; 

(iii)  generate  j,  a  random  variate  whose  probability 
vector  is  the  y_-th  row  of  P  =  Q/A  +  T; 

(iv)  generate  a  geometric  random  variate  with  param- 

eter  p(y  )  -  (I  -  q  i.e.,  P(C  »  n>  - 

n  y  n 

Cp(X  >)m_l0  ~  P(>’n)3*  an^  it  to  n; 

(v)  let  y  =  j ; 

(vi)  repeat  (iii)  through  (v)  until  n  ^  m. 

By  taking  r  replicates  we  can  estimate  PfY^  =  i)  using 
the  sample  proportion  >.(n)  of  replicates  for  which 
y^  =  i ,  i  t  S,  n  =  0,1,2, ...,m.  Using  Equation  (5), 
if  n  ^  T ( *  ,t). 


l  V«>  <«> 

n=0 

is  an  estimate  of  (t)  =  pfx(t)  =  i} . 

Algorithmic  Discrete-t ime  Method  (AD).  This  ap¬ 
proach  is  the  same  as  TD  except  that  the  discrete-time 

Markov  chain  (Y  ,  n  =  0,1,2,,..}  is  simulated  using 
n 

algorithms  to  compute  and  recompute  the  appropriate 
rates  and  probabilities  when  they  are  needed,  rather 
than  getting  then  from  precomputed  and  stored  values. 
The  probabilities  related  to  Y  are  then  transformed 
into  probabilities  related  to  the  continuous-time  pro¬ 
cess  X  by  using  the  randomization  formula.  Equation 
(8). 

The  two  algorithms  using  d iscrete-t ime  simulation 
and  then  randomization  are  attractive  because  they  of¬ 
fer  the  potential  for  variance  reduction  of  estimates 
based  on  the  simulation:  from  the  above  model  of  the 
cont inuous-t ime  Markov  process  it  is  clear  that  ran¬ 
domness  can  be  attributed  to  two  independent  sources: 

the  discrete-time  chain  (Y  ,  n  *  0,1,...  •  and  the 
n 

Poisson  process  !  X  ( t) ,  t  £  0-.  The  methods  TI)  and  AD 
simulate  Y  using  Monte  Carlo  and  treat  N(»)  determin¬ 
istic.  )]|v,  thus  there  should  be  less  variation  in  es¬ 
timators  based  on  1")  and  AD  than  in  the  sane  estima¬ 
tors  based  on  I'C  and  AC.  Crassmann  112|  has  also 
I'oi nt ed  this  out.  Hu-  question  is  whether  increased 
<omputat  ion  negates  the  benefit  of  the  decreased  vari¬ 
ant  e  of  t  In-  estimator. 

M!  I  !  I -l  tli!  I  .ON  RE  PA  !  HABIT.  ITEM  SYSTEMS 

"u 1 1  i -e* he  1  on  repairable  item  provisioning  sys¬ 
tem-.  »i*  vnri  1 1  i/.it  ions  of  the  classic  machine  repair 
o. ac  lonsidcr  a  svstem  consisting  of  two  bases 
md  i  -n  pot  .  I.idi  base  has  a  certain  number  of  nu- 

•  liiii'  -•  issi  ■■iii'il  l,'  it  ami  a  certain  desired  number  of 
it,,..,  vhi.h  should  be  I'perat  tug.  Machines  fail  (inde¬ 
pendent  I  v  .-t  i-.n-h  other)  alter  being  operated  for  an 

*  •  potu-nt  i  1 11 \  «1  i  st  r  i  but  t  d  length  of  time.  Then.-  arc 
r*  ”  i  *  r  shop-,  it  ,M,  h  has,  and  the  depot.  Vhrn  a 


machine  fails  it  has  a  certain  probability  of  going  to 
the  base  repair  shop  for  repair;  otherwise  it  goes  to 
the  depot  repair  shop.  Each  repair  shop  has  a  certain 
number  of  repair  channels.  Repair  times  are  exponen¬ 
tially  distributed.  If  there  are  more  machines  re¬ 
quiring  repair  than  repair  channels  at  a  given  repair 
shop,  a  queue  forms.  The  depot  repair  shop  stocks 
spare  machines  which  are  used  on  a  one-to-one  ordering 
basis  to  replace  failed  machines  coining  from  the  bases 
If  the  depot  spares  pool  is  empty  when  a  replacement 
is  required,  a  backorder  is  created  which  will  be 
filled  when  repair  is  completed  on  one  of  the  items  in 
depot  repair.  If  both  bases  are  awaiting  backorders, 
a  repaired  item  is  sent  to  the  base  with  the  maximum 
depot  backorders  (ties  are  broken  by  flipping  a  fair 
coin).  When  neither  base  is  awaiting  a  backorder,  a 
machine  completing  depot  repair  is  placed  in  the  depot 
spares  pool.  Thus  any  given  machine  may  be  in  any  of 
six  states:  failed  and  in  base  repair  shop  at  either 
base  (BR1 ,  BR2) ;  failed  and  in  depot  repair  shop  (DR); 
operational  and  at  either  base  (BUI,  BU2) ;  operational 
and  in  the  depot  spares  pool  (DU).  These  six  states 
and  the  possible  transitions  a  machine  can  make  be¬ 
tween  them  is  illustrated  In  Figure  1.  The  different 
parameters  of  the  system  are  described  in  Table  1. 


Base  1 


Figure  1:  General  schematic  for  a  two-base 
multi-echelon  repairable  item  system 


Table  1: 

Parameters  of  Multi-echelon  System 


Symbol  Definition 


BS.  Allocation  of  total  stock  to  Base  i  (operat¬ 
ing  machines  plus  spares),  1  =  1,2 

MS^  Desired  number  of  working  machines  at  Base  i 

BC.  Number  of  repair  channels  in  repair  shop  at 
1  Base  i 

FB.  Probability  machine  failing  at  Base  i  is  base 
repairable 

Aj  Failure  rate.  Base  i  items 
pj  Repair  rate,  base  i  items 

DS  Number  of  depot  spares 

DC  Number  of  depot  repair  channels 

Depot  repair  rate 

ij 


Steady-state  models  and  behavior  of  multi-echelon 
inventory  systems  are  presented  by  Sherbrooke  [13]  and 
Mucks tad t  [14],  A  method  for  computing  approximate 
transient  performance  measures  of  the  above  multi¬ 
echelon  system  is  presented  in  [15,16].  Gross  and 
Miller  [5]  compute  exact  transient  probabilities  using 
the  randomization  algorithm.  For  further  background 
on  nulti-echelon  inventory  systems,  see  the  above  ref¬ 
erences  . 

We  used  various  cases  of  the  above  multi-echelon 
system  as  tests  for  comparing  different  simulation 
methods.  Some  of  the  cases  considered  are  shown  in 
Table  2.  These  systems  are  initially  in  perfect  con¬ 
dition  with  no  failed  machines, 

COMPUTATIONAL  EXPERIENCE 

We  wrote  SIMSCRIPT  11,5  programs  to  simulate  the 
above  multi-echelon  inventory  model.  Five  programs 
were  based  on  the  five  Monte  Carlo  methodologies.  The 
event-scheduling  (ES)  program  is  written  in  the  style 
of  a  typical  SIMSCRIPT  program:  the  SIMSCRIPT  timing 
routine,  random  variate  generators,  and  future  event 
list  are  used.  The  bases  and  depot  are  modelled  as 
permanent  entities.  There  are  three  types  of  events: 
failure,  completion  of  repair  at  base,  and  completion 
of  repair  at  depot.  The  four  methods  for  Monte  Carlo 
simulation  of  Markov  processes  use  SIMSCRIPT  in  the 
style  of  a  general  purpose  language  rather  than  a  sim¬ 
ulation  language.  The  SIMSCRIPT  random  number  genera¬ 
tor  and  exponential  random  variate  generator  are  used. 
Geometric  random  variates  are  generated  by  a  routine 
that  either  uses  the  inverse  cdf  method  or  adds  the 
integer  part  of  an  exponential  random  variate  to  1, 
depending  on  the  value  of  the  geometric  parameter  (see 
[11]).  The  AD  and  TD  programs  use  £  *  .0025  to  set 
the  truncation  value  T(e,t)  in  Equation  (5).  The  low¬ 
er  tail  of  the  Poisson  distribution  of  N(t)  was  also 
truncated,  throwing  away  a  probability  of  at  most  .0025. 

The  deterministic  randomization  algorithm  (DET) 
was  programmed  in  FORTRAN.  It  is  described  by  Gross 
and  Miller  [5].  A  truncation  probability  of  c  =  .001 
is  used  in  Equation  (5)  for  this  algorithm;  this  guar¬ 
antees  that  all  computed  probabilities  have  an  error 
between  -.001  and  0.  Thus  we  are  implicitly  declaring 
an  indifference  to  this  size  of  error. 


Table  2: 

Parameter  Values  of  Mu  1 1 i -eche Ion  Cases 
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The  performance  measures  of  interest  are  availa¬ 
bilities  at  each  base  and  for  the  whole  system.  For 
i  =  1  and  2,  let  M^(t)  equal  the  number  of  operational 

machines  at  base  i;  if  M^(t)  *  MSj^  we  have  full  avail¬ 

ability  at  base  i.  Thus  we  define  the  following 
availability  measures: 

Avl(t)  =  PCM^C)  =  MS  ^3 

Av2(t)  =  PCM2(t)  -  MS23 

Avl2 (t)  =  p(m  (t)  =  MS1  and  M2(t)  =  MS23, 

for  t  ^  0.  We  shall  compare  the  different  simulation 
algorithms  by  the  efficiency  with  which  they  estimate 
Avl,  Av2,  and  Avl2. 

All  five  Monte  Carlo  programs  are  structured  sim¬ 
ilarly  for  sample  path  replication  and  data  collection. 
Estimates  of  Avl,  Av2,  and  Avl2  are  based  on  blocks  of 
100  replicates:  for  ES,  AC,  and  TC  100  independent 
replicates  are  generated,  the  availability  at  time  t 
of  each  noted,  and  the  sample  proportion  computed;  for 
AD  and  TD  100  independent  replicates  of  the  discrete¬ 
time  chain  (Y..,  n  ■  0, 1 , .  . .  ,T(t ,  t )  }  are  generated,  the 
availabilities  at  each  discrete-time  point  are  noted, 
the  sample  proportions  for  each  discrete-time  point 
computed  and  then  converted  to  estimates  of  availabil¬ 
ities  at  the  ront imimis-t ime  point  t  using  Equation 
(8).  To  get  more  accurate  estimates,  multiple  blocks 
(ot  100  r*pli«ates  each)  are  generated  and  the  esti¬ 
mates  from  each  are  averaged  to  obtain  single  esti¬ 
mates  tor  Avl,  Av2,  and  Avl2. 


I  he  »■  t  I  i>  icn.  v  rtu-usur 
m  t  o»  the  v.jr  i  an.  «  »t  t  he 
required  t  ♦  •  e%«-.  ute  the  pr 


measure  ot  a  program  is  the  prod- 
•  I  the  est  lm.tt or  and  the  CPU  time 
ttie  program: 


Hli-  lent  v  -  Variance,*  CPU  time. 

We  estimate  the  etfj.  iernv  bv  observing  the  CPU  time 
f  *  ►  r  ea*  h  exeiulitm  and  estimating  the  variance  (for  AD 
and  I'D)  or  using  the  exact  value  of  the  variance  (for 
KSf  AC,  ami  IC).  I'he  CPU  times  that  we  use  for  com¬ 
parison  Jo  not  itu  lude  the  time  required  to  compute 
the  generator  in  the  t able-look-up  method  (for  cases 
A,  B,  and  C  these  times  were  .82,  8.10,  and  31.24  sec¬ 
onds,  respectively)’  thus,  when  making  final  compari¬ 
sons  these  times  must  be  added  as  overhead.  The  DET 
algorithm  is  used  to  compute  the  exact  values  of  Avl, 

A v2,  and  Avl 2.  For  FS,  AC,  and  TC  the  variances  of 
estimators  based  on  100  replicates  are  Avl (1-Avl ) / 100, 
Av2 ( 1-Av2) / 1 00 ,  and  Av 1 2 ( 1 -Av 1 2 ) / 100,  respectively; 
thus,  these  values  are  computed  directly  (by  DET)  and 
the  only  purpose  of  the  Monte  Carlo  programs  is  to  ob¬ 
tain  estimates  of  CPU  execution  times:  10  blocks  (of 
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Comparison  of  Simulation  Methods 

We  expected  certain  patterns  in  the  performance 
comparisons  of  the  different  estimators.  Using  "table- 
look-up"  should  be  more  efficient  than  "algori thraic" 
if  the  time  to  compute  the  table  is  ignored.  AC  should 
be  slightly  better  than  ES  because  they  use  similar 
logic  but  AC  does  not  spend  time  managing  an  event  list. 
We  expect  the  "randomized  discrete-time"  approach  to  do 
better  for  large  t  than  small  t  because  of  T(c,t)  in 
Equation  (5):  the  rate  A  of  the  underlying  Poisson  pro¬ 
cess  is  computed  as 

A  *  MSt*X  +  BC1*u1  +  MS2*a2  +  BC2*U2  +  DC*pD. 

The  expected  number  of  events  of  N  (and  discrete-time 
values  for  Y)  in  [0,tj  is  At,  but  the  "randomized 
discrete-time"  approach  must  simulate  Y  for  0  ^  n  ^ 
T(e,t).  The  relative  extra  time  (j(e,t)  -  At}/At  can 
be  large,  especially  for  small  t.  Thus,  while  expect¬ 
ing  a  smaller  variance  using  AD  or  TD,  we  also  expect 
larger  CPU  times,  especially  for  small  At.  To  get  a 
better  idea  of  how  the  methods  perform  we  must  look  at 
some  comparison  data. 

Markov  Simulations  vs.  ES.  The  efficiencies  of 
the  four  Markov-oriented  Monte  Carlo  simulations  are 
Compared  with  ES  in  Table  3.  All  four  are  apparently 
superior  to  ES.  Note,  however,  that  AC  is  about  the 
same  as  ES  for  large  times.  The  others  (AD,  TC,  and 
TD)  require  significantly  more  programming  effort  than 
ES.  (The  ratios  in  Table  3  are  estimates  of  the  ratios 
of  efficiencies.  They  are  approximate  because  CPU 
times  are  approximate  and  because  the  variance  is  esti¬ 
mated  for  AD  and  TD.  Under  the  hypothesis  of  unit  ra¬ 
tios  and  deterministic  CPU  times,  these  two  sample  ra¬ 
tios  have  a  x?/df(60)  distribution.  Since 
XJ/df  (60;  .05)  =  .720  all  the  ratios  for  AD/ES  and 
TD/ES  are  statistically  significant.) 

Algorithmic  vs.  Table-look-up.  The  ratios  of  ef¬ 
ficiencies  AC/TC  and  AD/TD  are  shown  in  Table  3.  There 
is  a  clear  pattern  of  superiority  of  "table-look-up" 
over  "algorithmic."  Note,  however,  that  the  "table- 
look-up"  is  generally  more  difficult  to  program.  (Un¬ 
der  the  hypothesis  of  equal  efficiencies,  the  ratio 
AD/TD  has  an  F(60,60)  distribution.  Since  F(60,60; .95)= 
1.53,  most  of  these  values  are  significant.) 


Table  3: 

Comparisons  of  Markov  vs.  ES  and  A  vs.  T 


100  replicates  each)  were  computed  for  this  purpose. 

For  AD  and  TD,  60  blocks  of  100  replicates  each  were 

Kff  ic 

iency  Ratios  of 

Estimators  of 

Avl  2 

computed,  yielding  60  independent  estimates  of  the 
availabilities;  the  variance  of  these  60  estimates 

AC/ES 

AD/ES 

TC/ES 

TD/ES 

AC/TC 

AD/TD 

about  the  true  availability  (computed  by  DET)  was  com¬ 
puted  and  used  as  an  estimate  of  the  variance  of  the 

A 

3 

.58 

.72 

.33 

.39 

1.74 

1.83 

availability  estimator  for  each  of  the  three  avail- 

A 

14 

.87 

.56 

.50 

.45 

1.74 

1.25 

abi 1 i t ies . 

A 

95 

.95 

.52 

.55 

.26 

1  .74 

2.05 

We  initially  simulated  three  cases  of  the  multi- 

B 

I 

.  32 

.39 

.18 

.30 

1  .81 

1  .  30 

echelon  inventory  system  (Cases  A,  B,  and  C  in  Table 

2).  These  have  typical  parameter  values;  state  spaces 

B 

6 

.70 

.55 

.52 

.41 

1  .  34 

1  .  34 

of  375,  3960,  and  15075  states,  respect  Lvely ;  and  time 

B 

40 

.90 

.38 

.58 

.22 

1  .55 

1  .72 

points  which  encompass  transient  and  steady-state  be- 

.63 

.63 

.40 

.46 

havior.  We  then  ran  additional  cases  (including  D 

through  L)  to  gain  more  insight  into  the  behavior  of 
the  "randomized  discrete-time"  approaches  (AD  and  TD) . 

C 

30 

.90 

.55 

.54 

.21 

1.67 

2.68 

v;  /  v  y.v  -  vw-AVaV-V7  > 


Continuous- tine  vs.  Randomized  l)iscrete-t ime :. 

Table  4  presents  comparisons  of  variances,  CPU  times, 
and  efficiencies  for  AI)  vs.  AC  and  for  TD  vs.  TC.  We 
see  that  the  "discrete-time"  method  always  produces  a 
significant  decrease  in  the  variance  of  the  estimators 
of  availability,  but  also  a  large  increase  in  CPU  time. 
Resulting  efficiencies  are  sometimes  improved  (made 
smaller)  and  sometimes  not.  For  Cases  A,  B,  and  C  we 
see  a  pattern  whereby  the  "randomized  discrete-time" 
method  works  well  for  large  t  but  not  for  small  t. 
Initially  we  believed  this  pattern  was  caused  only  by 
the  value  of  A  =  £t(c , t)-At J/At  discussed  earlier. 

But  further  investigation  (Cases  D  through  I,)  showed 
that  while  a  large  value  of  A  usually  doomed  the  ran¬ 
domization  approach,  a  snail  value  of  A  was  no  guaran¬ 
tee  of  success.  For  Cases  A,  B,  and  C  the  randomiza¬ 
tion  approach  works  well  for  cases  where  the  process 
has  reached  steady-state  equilibrium  behavior  but  not 
during  initial  transient  periods.  It  is  interesting 
to  note  that  Hordijk,  lglehart,  and  Sehassberger  I  10] 
used  regenerative  steady-state  estimates  based  on  the 
uni  form! zed  embedded  chain  Y  and  achieved  approximate¬ 
ly  the  same  amount  ol  variance  reduction. 


Consider  the  behavior  of  base  1  in  Cases  .) ,  K, 
and  I.:  We  have  3,  5,  and  17  state  Birth  and  Death 
processes,  respectively.  They  are  all  in  steady  state 
and  there  is  no  activity  at  the  depot  or  base  2.  We 
get  approximately  the  same  improvement  using  randomiza¬ 
tion,  suggesting  that  the  size  of  the  state  space  may 
not  have  much  effect. 


Consider  the  behavior  of  base  1  in  Cases  F,  G,  and 
H:  In  each  case*  as  t  increases  the  availability  at 
base  1  remains  very  close  to  I. 00  and  then  decreases 
rapidly  to  close  to  0.0  (around  t  =  21  for  F  and  t  =  15 
for  G  and  H) .  Thus  we  are  estimating  Avl  at  nonequilib¬ 
rium.  Note  that  significant  improvement  in  efficiency 
occurs  in  Case  G  but  not  in  Cases  F  and  K,  and  further 
note  that  there  is  no  activity  in  the  rest  of  system  G 
and  much  activity  in  systems  F  and  H.  System  G  at  t  = 

15  is  the  only  example  we  have  found  where  the  random¬ 
ized  discrete-time  approach  improves  efficiency  in  a 
transient  situation;  this  example  is  characterized  by 
an  almost  deterministic  enbedded  Markov  chain  and  a 
moderate  value  of  £t(?.  , t)-At)//.t . 


Consider  the  behavior  at  base  2  in  Cases  D,  E,  F* 
and  H:  In  0  and  F  the  number  of  good  machines  is  a 
two-state  Birth  and  Death  process,  in  K  and  H  a  three- 
stale  Birth  and  Death  process.  Initial  transients  are 
very  strong  in  Cases  D  and  F.,  while  in  Cases  F  and  H 
the  systems  have  renehed  equilibrium  steady-state  be¬ 
havior.  In  Gases  F  and  l!  the  randomization  improves 
efficiency  hut  in  D  and  F  it  makes  it  worse,  in  spite 
ol  the  fad  that  variances  were  decreased. 


Comparisons  of  Cont i nuous-t ime  vs.  Discrete  Tine 


Ratios  for  Estimators  of  Performance 


Perf.  Variance  CPU  Time  Efficiency 
AD/AC  TD/TC  AD/AC  TD/TC  AD/AC  TD/TC 


ES  vs.  Deterministic  Randomization  Algorithm.  In 
order  to  compare  an  exact  method  (DF.T)  with  a  Monte 
Carlo  method  (ES)  we  modify  the  comparison  procedure. 
Consider  confidence  interval  estimation  of  p,  the  pa¬ 
rameter  of  a  Bernoulli  population:  In  general,  a  95 % 
confidence  interval  with  half-width  of  .03  requires  a 
sample  of  size  1068.  In  general,  a  99%  confidence  in¬ 
terval  with  half-width  of  .01  requires  16590.  In  order 
to  compare  ES  to  DET,  a  desired  confidence  interval  is 
specified  such  as  one  of  the  above  two  and  the  compari¬ 
son  of  ES  to  DET  equals  the  fraction  of  CPU .DET  which 
ES  needs  to  generate  the  necessary  sample  size.  Table 
5  gives  such  a  comparison.  It  can  be  seen  that  the 
comparison  depends  on  the  desired  level  of  accuracy; 
for  example,  to  estimate  Avl 2  for  system  C  within  i.03 
with  95%  confidence  would  require  at  most  only  19%  of 
the  CPU  time  required  by  DET;  but  to  have  ±.01  accuracy 
with  99%  confidence  may  require  triple  the  CPU  time. 

The  deterministic  algorithm  is  preferred  for  smaller 
state  spaces  or  when  a  high  degree  of  accuracy  is  need¬ 
ed.  llonte  Carlo  simulation  is  generally  infeasible  for 
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.049 

I  .  35 

1.56 

.057 

.076 

C  15075  5  ,059b 

332.30  5536.70 

.19 

2.97 

C 

15 

Av  1 

.076 

.092 

1  .49 

1  .86 

.113 

.171 

C  15075  30  .2651  1604.20  60*49.70 

.  18 

H 

15 

Avl 

.59 

.74 

1.33 

1.35 

.  787 

.998 

K 

15 

Av2 

.063 

.051 

1.33 

1.35 

.084 

.069 

* Average  time  per  re 

pl irate. 

J 

12 

Av  1 

.  n 

.15 

1.88 

2.78 

.249 

.415 

‘r  equals  the  number 

of  ES  replicates 

execut e 

d  i  n 

K 

12 

Avl 

.  1  7 

.21 

1.64 

i  .81 

.28  3 

.372 

tine  required  for  DKT 

execu  t ion . 

L 

16 

Avl 

.40 

.35 

1.21 

I  .21 

.483 

.428 

:: Ratios  of  C P U  t  i me s 

required  to  achic 

»ve  at  1 

east  95' 

12b 

confidence  interval  width  *.03  and  <19‘  confidence  in¬ 
terval  width  *.01,  respectively. 
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The  stiffness  of  the  system  plays  a  major  role  in 
the  comparison  of  ES  vs.  DET.  If  the  diagonal  of  the 
P,  the  transition  matrix  of  Y,  has  at  least  one  zero 
element  but  many  or  most  of  the  other  elements  are 
close  to  1,  Y  will  have  many  null  transitions.  Ran¬ 
domization  (DET)  loses  much  of  its  efficiency  for  such 
systems;  see  Miller  [7]  for  a  discussion  and  a  modi¬ 
fied  algorithm. 

CONCLUSIONS 

It  is  difficult  to  accurately  compare  different 
algorithms.  CPU  times  are  machine,  language  and  com¬ 
piler  dependent.  In  a  time-share/multiprogramming  en¬ 
vironment  these  timings  tend  to  exhibit  randomness. 
Behavior  of  the  algorithms  is  highly  dependent  on  the 
system  being  analyzed.  So  comparisons  such  as  these 
made  in  this  study  can  only  serve  as  rough  guidelines. 
In  view  of  this  we  can  still  draw  some  conclusions: 

Using  a  table- look-up  method  significantly  im¬ 
proves  efficiency,  but  it  involves  considerable  addi¬ 
tional  programming  and  cannot  be  used  with  systems 
which  have  huge  state  spaces.  The  randomized  discrete¬ 
time  approach  does  not  seem  to  improve  efficiency  ex¬ 
cept  when  the  system  has  reached  steady  state,  and  it 
requires  considerable  extra  programming.  The  deter¬ 
ministic  randomization  algorithm  works  well  for  sys¬ 
tems  that  have  a  small  to  moderate  number  of  states 
and  are  not  too  stiff  and  when  a  high  precision  is 
desired . 
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Decision  Tree  Models  Using  Monte-Carlo  Analysis 
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Decision  Tree  Models  lAinq  Monte-Carlo  Analysis 


Abstract  \ 

A  decision  tree  is  a  graphic  representation  of 
the  process  of  evaluating  the  value  of 
alternatives  leading  to  an  opromum  solution.  The 
process  of  decision  analysis  involves  (1)  graphic 
description,  (2)  assigning  valuiXand  probability, 
(3)  backward  induction  and  (4)  decision  making  or 
selection  of  the  path  through \the  tree.  A 

computer  is  effectively  used  to\  solve  the 

algorithms  of  Step  3,  the  backwani  induction. 
The  problem  of  reliability,  howeveV  is  not 

diminished  by  the  high  precision  of.  computer 
calculations.  The  results  can  be  no  belter  than 
the  best  guess  values  assigned  to  th\  various 
probabilities  and  outcomes.  Monte  Carlo  analysis 
provides  a  tool  which  can  utilize  the\  range 
between  the  best  case  and  the  worst  caseXwith 
best  guess  and  variability  to  quantify\the 
probability  of  achieving  anticipated  outcomes. 
The  added  dimension  of  probability,  significantly 
improves  the  reliability  of  the  decision  tree 
analysis.  A  methodology  is  proposed  for 
computer  model  to  solve  decision  tree  models 
using  Bonner  &  Moore's  Monte  Carlo  analysis, 
program,  PAUS.  A  standard  variable  naminc/ 
convention  is  offered  which  facilitates 
programming  decision  tree  models  and  a  samsfle 
PAUS  program  is  presented  as  a  guide  for  sol/ing 
decision  tree  models  using  Monte  Carlo  analysjfc. 

Overview  / 

Probabilistic  dynamic  programming  modads  sound 
significantly  more  important  than  decnsion  tree 
analysis  but  substantially  mean  the  srame  thing. 
People/managers  make  decisions  all  toe  time.  In 
response  to  the  morning  alarm  clock :/  (a)  get  up 
and  go  to  work,  (b)  call  in  sicyand  go  back  to 
sleep,  or  (c)  call  in  a  vacation  Aay  and  go  to 
the  beach.  The  logical  manner  makes  business 
decisions  based  upon  the  finanMal  value  of  the 
alternatives.  A  decision  Vree  analysis  is  a 
technique  of  quantifying  the  values  of  various 
possible  outcomes  and  /the  probability  of 
achieving  those  outcomes/ 

When  represented  in  /nathematical  notation,  a 
decision  tree  mode /  is  awe  inspiring.  However; 
if  represented  \x/  plain  business  language,  a 
decision  tree  \ys  simple  graphic  representation 
of  the  decision  Process.  The  computer  is  a  tool 
for  solving  nwthematical  relationships  that  have 
been  defined  isy  people.  A  manager  reaches  a 
decision  in/the  alarm  clock  problem  based  upon 
many  factor/  work  load  at  the  office,  other 
managers  Attendance,  personal  attitude,  the 
availability  of  sick  leave  or  vacation,  and  many 


Iempt  to  define  all  the 
te  to  the  value  of  each 
te  an  extremely  complex 
e  decisions  involving  many 
relative  ease  while  a 
Ive  the  same  problem  would 
.  Proper  model  design 
or  capability  of  the  human 
gical  work  and  let  the 
boring  repetitions  work 
machine).  The  decision 
llows  is  large  but  each 
erstand. 

rules  which  guide  the 
n  tree  solution.  (1)  A 
s  represented  by  a  square 
ts.  The  logical  manager 
deterministic  path  which 
turn.  (2)  A  probabilistic 
d  by  a  circle  in  the 
value  of  a  probabilistic 
sum  of  each  path's 
probability  multiplied  by  its  value.  (3)  The  sum 
of  the  probabilities  of  each  path  from  a 
probabilistic  decision  point  must  be  1.0. 

The  application  of  these  three  rules  makes  the 
solution  of  the  decision  tree  a  simple  (but 
tedious)  task. 

s.  Re  1  iabil  ity 

S\nce  the  arithmetic  calculations  involved  in  the 
baskward  induction  are  relatively  simple,  a 
computer  can  easily  be  programmed  to  compute  the 
results.  Unfortunately,  the  accuracy  of  computer 
arithmetic  infers  a  degree  of  reliability  which 
is  notNustified.  The  value  of  the  answers  can 
be  no  \  better  than  the  reliability  of  the 
estimates'^  of  the  costs,  revenues  and 
probabilities.  Since  the  computer  calculates  the 
backward  infection  values,  Monte  Carlo  analysis 
can  be  appYied  which  results  in  an  unlimited 
number  of  \solutions.  Using  Monte  Carlo 
techniques;  cqst,  revenue  and  probabilities  can 
be  defined  as  s\atistical  distributions  as  well 
as  discrete  Xoints.  The  results  increase 
reliability  by  \dding  the  dimension  of  the 
probability  of  achieving  the  various  decision 
point  values.  \ 

Monte  Carlo  Analysisy 

Bonner  and  Moore  have\  developed  a  Monte  Carlo 
modeling  tool  called  PMS.  The  PAUS  model  allows 
the  various  values  of  th\  decision  paths  to  be 
represented  as  any\  one  of  22  different 
distribution  types.  TheXmost  conrnon  financial 
distribution,  the  BETA  format  requires  the  user 
to  specify  the  minimum,  maxiYum  and  most  likely 
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